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SUMMARY 


A new explicit-implicit method for solving Navier-Stokes equations has been 
studied. The method has second-order accuracy in space and time, preserves the con- 
servation form, and is much less complex than other implicit methods since it 
requires no block or scalar tridiagonal inversions. It is used to solve a complex, 
two-dimensional, steady-state, supersonic-flow problem. A discussion is given of the 
computational efficiency of the new method and of the quality of the solution 
obtained from it at high Courant-Friedrich-Lewy (CFL) numbers. Modifications are 
discussed and certain observations are made about the method which may be helpful in 
using it successfully. 


INTRODUCTION 

Many numerical methods have been developed in the past several years to solve 
compressible Navier-Stokes equations, but the emphasis recently has been to develop 
implicit methods which are not subject to the conventional explicit stability condi- 
tions. However, the implicit methods are still restricted in time-step size by accu- 
racy and stability of the solution. Also, their programming complexity and computing 
time per time-step are much greater than those of the explicit methods. Recently, 
MacCormack (ref. 1) has developed a new method which is an implicit analog of his 
unsplit explicit method (ref. 2). This new method removes the explicit stability 
conditions, thus allowing larger time-step size to advance the solution. It can be 
added easily to existing codes which use the e^qjlicit method of reference 2. 

Further, the new method has second-order accuracy in space and time, preserves the 
conservation form, and requires no block or scalar tridiagonal inversions. Refer- 
ence 1 describes the method in detail and discusses the computational efficiencies 
obtained from it for a test problem. 

The gain in computational efficiency from the new method and its simplicity to 
be programmed make it very attractive over the other complex implicit methods. This 
method has been incorporated into an existing code developed in reference 3 which 
uses the unsplit explicit method of reference 2. The code is operational on the 
Control Data CYBER 203 vector processing computer. The purpose of the present report 
is to indicate modifications that were made and to make certain observations which 
may be helpful in using the new method successfully. The method has been applied to 
a two-dimensional supersonic-flow problem involving shock and expansion waves and 
their interactions with each other and with the boundary layer. A discussion of the 
computational efficiency and the quality of the solution with increasing time-step 
size is also presented. 


SYMBOLS 

CFL Courant-Friedrich-Lewy limit 

c velocity of sound 


specific heat at constant vol\ime 



e 

h 

M 

Npr 

n 

P 

R 

T 

t 

At 

u 

V 

x.y 

Ax 

Ay 

Y 

P 

Subscripts: 


total energy per unit volume 

static enthalpy 

Mach number 

Prandtl number 

number of time-steps 

pressure 

gas constant 

temperature 

time 

time- step size 
velocity in x-direction 
velocity in y-direction 
Cartesian coordinates 
grid spacing in x-direction 
grid spacing in y-direction 
ratio of specific heats 
effective viscosity, \x^ + \x^ 
density 


Z laminar 

t turbulent 

" free stream 


GOVERNING EQUATIONS 

The two-dimensional Navier-Stokes equations in conservation form are used to 
describe the flow field* These can be written as 
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In order to complete the set of governing equations, the equation of state 
p = pRT is used. In the aforementioned equations, \x is the sum of laminar viscos- 
ity and turbulent viscosity. The laminar viscosity for air is calculated from 
Sutherland's formula. The turbulent viscosity is calculated from an algebraic, two- 
layer, eddy-viscosity model due to Baldwin and Lomax (ref. 4). 
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METHOD OF SOLUTION 


The governing equations are solved by a new method, developed MacCormack 
(ref. 1), which is an implicit analog of his unsplit e:q>licit method (ref. 2). The 
method consists of a predictor step and a corrector step which can be written as 


follows (nomenclature 

used is 

same as 

that 

in ref. 1 ) : 

Predictor step: 
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Corrector step: 
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Each step contains two stages. The first stage uses the explicit method which is 
subject to restrictive explicit stability conditions. The second stage removes these 
stability conditions by transforming numerically the equations of the first stage 
into an implicit form. It is seen from the equations that the method requires the 
solution of upper or lower block bidiagonal equations. 

* 

Here/ I is the unit matrix and oU/ 6u , and so forth, represent the change 
in U with time— step; A = 9 f/3u and B = 9 g/9u are the Jacobians of F and G« 
The matrices |a| and |b| are the matrices with positive eigenvalues and are 
related to the Jacobians A and B. Their definitions are given in reference 1 and 
are not repeated here. This method is more efficient than existing methods for 
solving the equations of compressible viscous flow because, for regions of the flow 
in which At satisfies the explicit stability conditions, the method reduces to the 
simple explicit method, and, for other regions, block bidiagonal equations need only 
be solved rather than the block tridiagonal equations of the existing implicit 
methods • 
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BOUNDARY CONDITIONS 


The flow variables at the inflow boundary are held fixed at given free-stream 
values, whereas first-order extrapolation is used to obtain the flow variables at the 
outflow boundary. Reference 1 uses reflection-type boundary conditions which are 
not, in general, convenient to use. In the present calculations, no-slip boundary 
conditions are used along the solid surfaces. 

To start the implicit sweep in the i-direction (i varying from 1 to I), 

|a|6Uj . is required in the predictor step and |a|6u^ . is required in the 
corrector step. Both of these end fluxes in the i-direction are set equal to zero 
for the present problem because, at the inflow boundary, the conditions are held 
fixed and, at the outflow boundary, the mesh spacing is large enough so that matrix 
|a| vanishes. The implicit sweep in the j-direction (j varying from 1 to J) 
requires evaluation of either |b|6u^ j or |b|6u^ • These are obtained by 

calculating |b| and 6u at the boundary points. ^Matrix |b| is calculated by 
using the boundary values of the flow variables from the previous step, whereas 
6u at the boundary mesh points is obtained based on the changes in U at the adja- 
cent mesh points, again from the previous step. As an example, for no-slip and 
adiabatic-wall boundary conditions with 


Pi,1 " Pi, 2 

T. - = T, ,, 

1,1 i / 2 

the ( 6u ) . - can be written as 

X/ 1 

Ni.1 ■ (S)i,2 
Ni., - ■ 0 

Wi,i ■ - ’) 

where U^, U 2 , U^, and are the components of U. 


ARTIFICIAL DAMPING 

To maintain a stable solution, it was found necessary to use artificial damping 
in both the explicit and implicit stages of the method. The fourth-order damping, 
already present in the explicit code of reference 3, is retained for the explicit 
stage of the method. For the implicit stage, the damping term given in reference 1 
is used. This term is written as 

|6p/c^ - 6p| 

(At/Ay) [(Y - 1)/Y]P 


The contribution from this term becomes small as the steady state is approached. The 
method developed instabilities if either of the two damping terms was eliminated. 
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DISCUSSION OF RESULTS 


The modified code is applied to a two-dimensional test problem shown in fig- 
ure 1, in which the turbulent supersonic flow is calculated between two walls for air 
under perfect gas assumption. The upper wall is kept parallel to the free stream, 
whereas at the lower wall, the flow undergoes a 10® compression at 2 cm and then a 
10® expansion at 4 cm. The total length of the flow domain is 10 cm with an initial 
height of 2 cm. This problem is chosen since it involves most of the features of a 
complex supersonic flow, such as shock— expansion-wave interaction and shock— 
boundary-layer interaction which results in a separated region. Solutions are 
obtained at CFL values of 1, 15, and 22.5 for inflow conditions of M = 5, 

PcD =0.1013 MPa, and T^ = 293 K. The Reynolds number based on the free-stream 

conditions and the length of the flow domain is 11 10®. An assessment is made of 

the efficiency of the new method on the Control Data CYBER 203 vector processing 
computer and of the quality of the solution obtained from the method. 

The physical domain, shown in figure 1, is numerically transformed to a computa- 
tional domain. Since, for the present problem, matrix |a| vanished at all the grid 
points and the grid is not skewed in the x-direction, the numerical transformation 
did not affect the implicit steps. In general, matrices |a| and |b| need to be 
defined properly in the transformed coordinate system. A grid size of 51 x 51 is 
used in the calculations. To resolve the turbulent boundary layer, a stretching 
function is used in the transformation that allows concentration of grid points near 

the walls in the physical domain. The grid in the computational domain still remains 

equally spaced. The coefficient in the stretching function is chosen so that the 
first grid point away from the walls is located at approximately 4.5 |im. The 
computer code is written for the Control Data CYBER 203 vector processing computer. 
The explicit steps of the code are fully vectorized, whereas only a portion of the 
implicit steps could be vectorized. 

To start the solution, free-stream conditions are assumed initially at all the 
grid points except at the boundaries where appropriate boundary conditons are used. 
This starting solution works well for CFL values of 1 and 15, but for CFL = 22.5, 
negative temperature developed near the expansion shoulder in the first few time- 
steps. For this case, the solution is advanced with CFL = 15 for the first 
100 time-steps to establish a reasonable initial solution, and then the CFL value 
is increased to 22.5 over the next 200 time— steps. The solution is advanced in time 
until it reaches the time required by the free-stream to traverse the flow domain 
three times. As mentioned in reference 1, At/p(Ay)^ must be less than 0.5 to 
avoid any possible dependence of the steady-state solution on At. In the present 
calculations, this quantity exceeded 0.5 for CFL values of 15 and 22.5. For these 
cases, the time-step was successively reduced near the end of the calculations. 

It is found that the present problem at high CFL required the reversal of the 
order of differencing for the predictor and corrector steps from one time-step to the 

time— step? that is, if forward and backward differences are used in time— step 
n, then backward and forward differences should be used in time-step n + 1. Without 
this symmetric operation of the differencing, the solution at a CFL value of 10 
in the 27th time— step. Even at a low CFL value of 2, the quality of the 
solution is much better with the symmetric operation. 

In order to assess the quality of the solution with increasing CFL, pressure 
distributions on the upper and lower walls are compared at CFL values of 1, 15, and 
22.5. Also compared are the velocity and pressure profiles at three locations on the 
upper and lower walls. These locations lie ahead of, aft of, and in the separation 
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region. Figure 2 shows the pressure distributions on the upper and lower walls. It 
is seen that the lower wall pressures are almost identical at the three CFL values, 
whereas on the upper wall there are small differences near the 0.07-m location. This 
is the region where the shock from the lower wall separates the boundary layer on the 
upper wall. 

Figures 3 and 4 show the velocity and pressure profiles, respectively, at 
X = 0.05 m, which lies ahead of the separation region. It is seen that there are 
negligible differences in the profiles at various CFL values. Figures 5 and 6 show 
the velocity and pressure profiles, respectively, at x = 0.074 m, which lies in the 
separated region. The velocity and the pressure profiles at the lower wall again 
compare very well, whereas the velocity profile on the upper wall shows significant 
differences close to the wall at various CFL values. It is seen that the extent of 
the separation region is reduced at higher CFL values. The differences are also 
present in the pressure profiles on the upper wall. The separation region may have 
been affected by the artificial damping in the explicit and the implicit steps which 
is required to keep the solution stable. 

Figures 7 and 8 show the velocity and pressure profiles, respectively, at 
X = 0.088 m, which lies downstream of the separation region. Here again, the 
profiles compare well at various CFL values. 

The increase in computing time due to the addition of implicit steps is approxi- 
mately 50 percent on the Control Data CYBER 203 vector processing computer. It takes 
approximately 2 10" sec per grid point per time-step for the explicit method. 

This time increases to about 3.0 10 ^ sec with the addition of the implicit steps 

which could be vectorized only partially. Even with the increased computing time per 
time-step, there is a significant savings in overall computing time. As an example, 
the total computing time required is approximately 10 times less at CFL = 15 than 
at CFL = 1. 

As with the other implicit methods, the gain in the computing efficiency is 
expected to be higher with a highly stretched grid in one coordinate direction. 
Although reference 1 describes the method as unconditionally stable, it is not found 
to be so in the present calculations. For the present problem, solutions could not 
be obtained for a CFL value of 30 or above. 


CONCLUDING REMARKS 

A new explicit-implicit method has been applied to a complex, two-dimensional, 
steady-state, supersonic-flow problem by using Navier-Stokes equations. The method 
originally used reflection-type boundary conditions, but in the present calculations, 
no-slip boundary conditions are used. It is found for the present problem that the 
method requires the reversal of the order of differencing in the predictor and 
corrector steps from one time-step to the next time-step. It is also necessary to 
use artificial damping in both the explicit and implicit stages of the method to 
maintain a stable solution. A detailed comparison of the results at Courant- 
Friedrich-Lewy (CFL) values of 1, 15, and 22.5 shows that the flow field is predicted 
very well at high CFL values except in the immediate neighborhood of separation. 

For the present problem, solutions could not be obtained for a CFL value of 30 or 
above. 


The code is operational on the Control Data CYBER 203 computer. There is 
approximately a 50 percent increase in computing time per time-step due to the 
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addition of implicit steps in the explicit method, but the overall saving in the 
computing time is very significant. As an example, for a CFL value of 15, the 
computing time is reduced by a factor of 10 over a CFL value of 1 . It can be con- 
cluded from this study that the method has great potential in computing complicated 
fluid-dynamics problems since it is computationally efficient and is relatively much 
simpler than other implicit methods. Existing explicit codes using the unsplit 
MacCormack method can easily be modified for the new method. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
September 23, 1981 
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Figure 2.- Surface-pressure distributions. 
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Figure 6.- Pressure profile at x = 0.074 m. 
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Figure 7.- Velocity profiles at x = 0.088 m. 
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